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The purpose of this thesis is to develop an algorithm for segmenting images 
Semiotedwoy a men level Of noise wit diilerent characteristics. [n particular the 
images considered are composed of several regions describing different objects aid 
background. The algorithm described is based on a Markov Random Field (MRE) 
Medel Ol the image and uses Kalman Filterin@ (KF) techniques and Dvynanuc 
Progranuming (DP) 1n order to smooth within the regions. The theoretical background 
tor one dimensional and two dimensional data Which have different characteristics and 
simulation results are presented, with examples of synthetic data and underwater 


ina ces. 
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I. INTRODUCTION 


The objective of segmentation is to divide a given image ito meanimeglu regions 
or units, In many vision problems the first task is to distinguish objects to Jiscrutunate 
between Various obstacles. 

A general vision system for Artificial Intelligence applicuuions can be decomipesed 
in the following blocks: 


lL. Data acquisition: The picture of interest can be taken by suitable cameras for 
these purposes and can be digitized using proper software and hardware. 


tY 


Segnientation. Removing noise and dividing the image into significant regions. 


3. Image understanding: JJeciding which objects are present and classifving them 
ACCOR IMCHIO sIZe, spe secre, 


The general approach at the basis of image segmentation its based on the 
identification of different properties characterizing the regions of imterest. For 
example. several objects can be identified by their average intensity levels or by the 
textures on the surfaces. The task of an tunage segmentation algorithm 1s therefore to 
identify the regions from the vision signals which contain noise (due to the electronic 
equipment and other environmental factors such as turbid waters in underwater 
environments) or uninteresting details. For example, if it is desired to recognize the 
presence of a house in the picture, we need to segment the image by disregarding the 
GGors, WintloWs, cricks on the wall, etc: 

Although all the methods in the hterature for segmentation are well suited to 
niost of the cases of interest, they have poor performaiices i the presciice of aman 
level of noise. In the cases of noisy data statistical techniques are more suitable. In 
particular these are based on classical techniques of estimation, such as \faximum 
Likelihood (ML) or Maximum a Posteriori (MAP). Statistical estimation techniques 
rely on the use of distinct statistical models for the original image and for the 
disturbances. 

In this thesis we present a segmentation algorithin for images affected by a high 
level of noise, based on statistical techniques. In particular we use Markov Random 
Fields (MRF) as a model of the original image, and we assume White Gaussian noise 
as a model for the disturbance. The reason for the choice of SIRF is based on several 
considerations: 


1. They extend one dimensional well-known models to the multidimensional case. 


tJ 


They relate global statistics to local dependencies of the data. 


fy2 


They can model compact, distinct regions bv a suitable choice of parameters. 
The smoothing 1s obtained by a combination of MRF (for the transition between 
regions), Kalman Filtering Techniques (to filter within the regions) and Dynamic 
Progranimung (to determine the best sequence of edges which maxinuzes the likelihood 
function). 

Related to the works of segmentation, the algorithms of Derin, Elttott and Cristi 
pec (rcttiaty |Ref. 3], Besag fiefs) have shown the effectiveness of these 
techniques segmenting images corrupted by a high level of noise. Paratlel to these 
works, a combination of Autoregressive and MRF models has been presented by 
Therrien [Ref. 5] to segment textured images of terrain data. 

In the next chapter, the properties of modeling the original scene by Markov 
Random Field and choice of the parameters in Markov Random Fields are presented. 
Estimation algorithms for one and two dimensional data are given in Chapter IH. 
Finally, simulation results for different characteristics of data are the subjects of 


Giaprer |v. 


ff. STATEMENT OF THE PROBLEM AND MODELING USING SIRF 


The problem to be addressed in this chapter 1s low tormoedel the omoiiaie ean. 
and its estunation from the given (observed) data. As stated in the Introduction the 
model for the original scene is assumed to be a Markov Random Field (MIRF> whose 
properties are discussed in the next section. As, explained below siliceestame me 


based on a Bayesian approach and a Maximum a Posteriori (MAP) techniques. 


A. PROBLEM STATEMENT 


Let the original scene X = {X, ,{ be a random field on a finite two dimensional 


° 2, A : = . ° 
lattice L = {(Kytpe Z 2 tl Sk S N,.1 St S&S X,} where Zs tne scr or tes a 
assunung discrete values F,, F,. F3..... which are constant in regions R,. R,. R3.... with 


R, © Land R, n R, = 0 for i+). For example this might be the case of several objects 
MithechTerenitwmltensit.-1e.els- 
Also let X be corrupted bv an additive noise and modeled by a random field 
W= eae (k,t) € L. Furthermore, W will be assumed to be a White Gaussian@itela 
with zero mean and known variance o7. W is identically independently distributed 
(i.1.d.) and independent of X. 
Therefore. the observed image can be given bv a random field Y = eens aS 
Y (2.1) 


X W 


ec Ci See ie 


where (k,t) € L and a (.,.} 18 an arbitrary function. In our case we assume acdiine 
noise, and therefore the observed signal can be expressed as 


K, k,t 


Now, the problem is to estimate X from the data Y. Since we assume that only 
noisy observations are available, the estimate can be obtamed by maxunization of the 
a posteriori distribution of the scene with respect to x by a Bavésian opprogcm 


Therefore 
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where X is the estimate of original seene and Prey (sl) SLO tse t= fC), 
As it can be seen from Equation 2.5, the problem ts to maximize the ua postericr: 
probability of X given v. This form of Basesian estimation is known as Masui o 
Posteriori or MEAP estimation. 
Bayes factorization and the convenience of logarithmic operation Vields 
A A | | 

InP yin (¥ix) ae (tt OS LAX Py (sls) cegliilex 0) oe 
The two likelihood terms on the left and right hand sides of Equation 2.4 ure 
ieommuncuse. the amodelio: \ and of tiie distirbunce. The model! for the scene X Is 
assumed to be a Markov Random Field (IRF) which is formally delined by Besag 
{Ref. 3: p. 724] and by Elliot. Derin, Cristi and Geman [Refs. 6.7]. 


B. MARKOV RANDOM FIELDS (IRF) 
Pomorie bet Le be aime lattice Lo = {(kt): 1S k Ss Vote a aod 
nN, , defined as a subset of L so that (k,t) € n, , where (kt) € L and (1j) en, if and 


Caeeuetimt eG i. Then arandomiielid \ = {N, .} with the property 
2 hij t kw p p a 


PiX 
PIN 





AigaXi pe Ud) € Li = eo 
oor - ‘ 
Ni Nie (1.J) € Mey 


kt hat 





ket 


isa Markov Random Field with neighborhood n, .. 
Poa nov Random Picld can be illustrated as in Figure 2.1, where the statistics 
Oiieerc ment (kaminciedied bY a ~ depends on ifs newehbors indicuted by a “QO” 


only. Two simple neighborhoods are shown in Figure 2.2. These are: 


le retin): 0 < (i-1)° a (j-n) = |} fo 
and 
TH = {(1,n):0 < (i-1)* + (j-ny? = 2} ea) 


Relatively simple neighborhood systems as in Figure 2.2 are adequate mn modeling most 


scenes of tnterest. 





Ficure 2.) Dependence of Starkey Randomeniem 
P("|Everytiung else). 











(a) (b) 





Figure 2.2. Neighborhoods n, ‘i (a) and n, & (b). 


Although a MRF can be considered as a muludimensional extension of a 
Markov Chain, a difficulty exists in definng MRF. In fact the main concept at the 


basis of Markov Chains is the concept of transition probabiliues. These are subjects to 


nuld conditions and thev rely on the fact that an ordering can be deternuned in one 
dimensional problems. This is not the case in multidimensional problems where an 
ordering of points in general does not exist. Which presents a difficulty in the definition 
of transition probabilities. Because of this we have to approach MRFs in a different 
ameonmecie = larko. Chains. ihermvere, to deternune the structure for the joint 
PremavilliwecmmvtrRd models, we need some new definitions and theorens. In 
particular the joint probability is based on the concept of cligue given below: 

Definition; Given a Markov Random Field (MERE) with neighbors 1) a cljue is a 
Sema pixciommmich are neighoors of cach other (Rei. 6: p.19S]. The various types of 
cliques for is and ne Meson mierieier.>..) Based On tie deiimition of clijues. 
the Joint probability of a Markov Random Field with neighborhood ny can be expressed 
Dwene [Ollowine theorem: 

Theorem (Hanimersley and Clifford) [Ref 6: pp. 192-199]. Suppose X is an MRF 
such that Py(x) > 0 for all x, with neighborhood yn. Then Py(x) can be defined us 


following: 


l ie 
Px(eseh ®) Ces 
Where 
U(xy= VV (x) (2.9) 
ces 


5 is the set of all cliques as shown in Figure 2.3 and usually U(x) is defined as the 
energy function, V(x) the potential associated with clique c, and the partition function 
Z is defined as 


‘ 


which ts a normalizing constant that causes Py to be a consistent probability measure. 
On the basis of the above theorem we can arbitrarily assign a joint probability of 


MRFs by the potential functions V(x). The only constraint on V(x) is that it has to 





(a) 

















(Db) 





Figure 2.3 Cliques for Neighborhood Systems ni (a) and ns (b). 


be finite for alleweelt tacans tare ies onnamon V (x) is fixed by the structure of the 
lattice 2iven by Kimdemman, | Rew 
in order to understand the above result and definitions, the joint probability of 
anv Markov Random Fteld with neighborhood Nis can be written as 
InPy(x)= ie ar Be eae ae Sal InZ (27 
k,t kt kt 


(OF Ine Casein Pigtinesen a). 


Definition; A MIRE which has all nonzero probabihties as in [qguattun 2.11, 
namely, Py(x) > 0, is called a Gibbs Field (GT). The Gibbs medet can Fe used to 
model the spatial interactions between regions, in which a pixel and its neighbors lave 
sinular statistical properties. [In other words. tt models spatial continuity which 
characterizes the images we want to segment. Note that the structuic of P (Xx) ts 
relatively siniple for Thre (See ia OMmeelmelis CONDIENILY tiicreases Considerably tor 
larger neighborhoods. For that reason. in general onlv the lower sets Ns and eee are 
considered practically applicable, and even in Me Chdicce Vath conly One Or «(wo 
elements are taken into account. In addition, the Gibbs distribution ought aiso be 
used to model textures. In our problem we just want to model the fact that the most 
likely scenes have regions which are clusters of pixels, and the Markov Random Fieid 
minodel is adequate for this purpose. provided by¥ properly assigning the potential 
funetions. A particular model we are going to consider in this thesis is the M&F with 


neighbors ny! as in Figure 2.3 (a) and with potential functions defined as: 
, =— fy as 
cae ° oe 


Bifx, =x 

aia _ : = ek k t-] 4 
ew, ) = VAN, Ky eae _—. 
for any value of k and t and B a positive constant. The larger the parameter PH. che 
niore the joint distribution Py is peaked around smooth realizations. By this definiuon 


we can Write the joint probabilities as 


TN Pe 5 i) (Ky Xe} InZ (2.14) 
Kt 


where 


Pix = Xe 


. — >) ~ 
Q(X os Xeepe! = (23ley 


-1 otherwise 


Also, by the assumption of Gaussian noise with zero mean and standard deviation o, 


Equation 2.14 with Equation 2.4 yields the likelihood as 


Hela 


FAR, Xe es 


| 
x (SIX) + InPy (x)= = Sr “= Bip 
Fike * 


InPy- 





Xk ot kel ot 


Kp in noise, the noise model W(®, 6) and the Gibbs field 


model Equatron 2.14 with parameter B, the estimate of the original scene is given by 


Hence, given the data y= {v 


os h th 
x= cree) such tnat 
are =X ee Oe (X Vee Nee )} Pale 
Igo KA Kt iy ke tel Sk tke Lt a 
ket kt 


Is minimial. 

In particular there ts no need to undertake the difficulty of calculating the 
partition function Z since it is a normalizing constant. So, it was omitted in Equation 
2.17. Approaches to nunimization of Equation 2.17 can be based on relaxation 
[Ref. 3: pp. 730-732] or deterministic [Ref 4: pp. 266-267] or Dynamic Programming 
(Ref. 2: pp. 44-45 and Ref. 8: pp. 12-13]. In the next chapter a different approach to 


nunimuzation of Equation 2.17 based on Kalman Filtering techniques will be presented. 


C. CHOICE OF THE PARAMETER IN THE MRF 

The difficulty in using the Grbbs Distribution as a region model is to estimate the 
parameters of the model from specific realizations. Some methods have been used 
previously to estimate the model parameters. For example, Besag [Ref. 6: pp. 211-212] 
suggested the coding method where the parameters are determined from subsets of data. 
This requires solution of a set of nonlinear equations. Another example of the 
paramieter estimation method its given bv Derin and Elliot [Ref. 2: pp. 43-44] which 
uses standard ltnear, least-squares estimation. This has been proved to be effective mn 
modeling textures by GFs. A different approach was used by Geman [Ref. 3: pp. 
727-729] which deftned a simulated annihiling, where stochastic relaxation was 
combined with monotonic increasing parameters. 

For this thesis, the parameter B of the model in Equation 2.17 was set bv trial 
and error until a reasonable ftitering was reached. The optimum value of the 
parameter [ is smaller for high values of stgnal to noise ratio. So, it is necessary to 


modify this parameter as the signal to notse ratio changes. 


(2.10) 


I. SEGMENTATION AND S\IOO THING OF ONE DINIENSIONAL AND 
TWO DIMENSIONAL SIGNALS 


The content of this chapter is the seginentation algorithms for one dimensional 
and two dimensional signals. Following a MMaxinuen a Posteriori approach, esumate of 
the onginal data is obtamed bv mininuzation of the cost funetion given by Equation 
ieee coe provide this, a stmoothine Ttechimigue is devised which 1s based on a 
combination of Kalman Filtering (KF) and Dvnamic Programnung (DP). Dynanuc 
Paeciamunmimneiseused=to deterimime the best seuuence of edges which maxniuzes the 
hkelihood function and a detailed algorithm is presented in Appendix A. Kalman 
itemise ised to filtereuion the tegions. This filtering technique is preterred 
because the Kalman filtering is the best linear optimal estimator. I[f the noise in the 
data 1s assumed to be Gaussian, the Kalman Filter gives the minimum Variance 
estimate of the original scene. In particular tt evaluates the conditional mean of the 
original data given the past observations (measurements). If the assumption of 
Gaussian noise is removed, the KF vields the best linear estimate. 

In the next subsection. the segmentation algorithm for a binary sequence (1e., 
one which assumes onlv two levels) is presented and the result will be extended to 


general multilevel cases. Related computer programs are included in Appendix B. 


A. BINARY SEGMENTATION USING DYNAMIC PROGRAMMING 
ALGORITHMI 


Suppose that it is desired to segment a binary sequence having logic levels “1” 
and “O” which are related to two intensity levels (grav levels) corrupted by an additive 
noise. The noise is considered as a zero mean Gaussian noise with a known standard 
deviation 6. In the general image processing extension, we can consider a logic “O” as 
a background and a logic “1” as an object. “An example of the original sequence and 
the noisy sequence is given tn the next chapter. 

Let x, be the logic values of the original data of intensity levels, which «ssigns 
real numbers F(0), F(1) to the regions denoted by logical zeroes and ones respectively. 


Then the norsv data can be given as 
eae) > W, (Sale) 


where W, ~ N (0,6). We can define the signal to noise ratio (SNR) as 


we 


+N 


S\R = ————— ee) 


\ 


[n this section we assuine to know the levels [(0)) fl)” andethe noise Varian ccm 
Here, the assumption of F and 6 known ts not restrictive. since various algorithms 
exist in the literature which enable one to estimate means and variance of nuxtures of 
Gaussian populations. We have seen in Chapter 2 that a MIRF can be arbitrarily 
assigned by the potential function V. In particular we can model spatial continuity by 
assigning hich probabilities to smooth signals. This can be achieved by assigning the 


likelihood function in the following form: 





Sy ~, 
(kee a : Wy.-F(x.)|- ee 
r=O0 1=0 
where 
lif x,=x 
! : ee 


g(X 4X4) = : 
‘ -! otherwise 


The parameter Bh models the sinoothness of the original data. In recursive form 


the likelihood function can be written as: 
| I 5 t 
Ee (Xp eXpreeeXy $= AKG Xp eeXy) + BECK, 4 PX) 5 or TK fp US sell (3.9) 


with €, = 0. The estimate x is obtained by maximizing € over all possible realizations 
x. Applying Dynamic Programming(DP) algorithm in Appendix A, the maximum of 


Si, is obtained. Results of segmentation of noisy binary iunage ts given in Chapter 4. 


Ls 


Cee OUIING SOF PIECEWVMSE CONSTANT SIGNALS IN ONE 
DINIENSION USING KALMAN FILTERING 


The algorithm presented in the previous section for binary data is suntabe for a 
siall tuumber of levels. Although it can be generalized to anv number of levels. its 
complexity increases considerably and becomes unfeasable in multilevel situations. So, 
for one aiid two dimensional observations which have more than two intensity levels, 
Peewmeccmo Searcl) for a dilierent approach, 

Psementioned im the first Section of tls chapter. a new algorithin wlich uses the 
Kalman Filtering techniques and Dynanuc Progranuning is developed. In appheations 
filtering 1s concerned with the extraction of signals from noise. [f the observation and 
original scene are modeled by linear stochastic models, a solution to the general 
filtering problem can be obtained bv use of the Kalman Filter. 

I Wem on tile scscUNiedmmiccewise Gaaticteristics Oo: the data, the realizanons of Y 


can be modeled by the state-space models 
ae n) 


t 


where w, is Gaussian zero mean itd. with standard deviation 6, and vy, 1s nonzero at 
the edges between regions onlv. By defining the ornginal data X mn one dimension as 
eee Le (ioe. the one duncnsional Markov model corresponding to 


MRE has joint density 


InPy(x)= BY 9(x,,.3,.1) - InZ (3.59) 
teL 


and the likelihood function to be maxinuzed like Equation 2.16 1s 


I : 
{(x)= - Ser” + BY o(x,.x,_,) (So) 
teL teL 


Where the first term on the right hand side of this equation comes from the notse and 


the second term 1s the potential defined before. 
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The Kalman Filter for the state-space model in Equation 3.6 and 3.7 1s defined 


bv the recursion Equation 3.11 as 
a ‘ Va Oe “ LG: 
Soy <> See (9.10) 


Where K, is the Kalman gain. When an edge is detected tie’ cain nesdss temo 


reinitialized at that point. To do this, define a binary sequence 
90 bates (Sal iy 


which has a logic “1” at the edges between the regions and a logic “0” in the regions. 
Just to give an exaniple, let’s say T, = [t,. t, - l] and T, = [t, , t, - 1] where T,. T, are 
adjacent regions and t, < t, <t,. We can define ¥,, = V,. = V,, = 1 as shown in 
FiQire- sei linesemenan 


t 


{ lift- eT, andté T, for some k 
—— 


QO otherwise 


| 
| 
Re 
| " : 
| , : 
Viq =! | 
‘ | 


Then the estimated data can be determined from the noisy sequence wing Dyaauion 
Belvrand 3.11 as 


eee Niles, el 
where 
; l ‘ 
x, =— ie 
. 
t 
acl: 
in =! ; 
i (ors) 


t., + | otherwise 


Filtering Equation 3.13 vields the estimate of x as 


1 
l t 
= ae yy 5 sar: 
i : oj SSIS 


t j=l 


Where t, stands for the distance of t from the nearest left edge. Equation 3.16 
expresses the Kalman Filter as an averaging of the measurements. 
In the case that the edges are not known, observe the following: 


: : : c LS 6 i = 
|. Given the noisy data v for anv binary sequence m, the estimate X is weil defined 
‘ = . : A 
Poeequationvaaks = 5.15, call this estiinate x(m); 


2, The g(.,.) term in the Markov model Equation 3.8 and the likelthood function 
Eee Jecaiebeiiterpreted aS a penalty to the edges. From the 
dependence of x on the sequence m stated above it can be written. 


mi = 0 


-lifm, = 1 


Therefore we can define the function 


| {ifm = 0 
NU) = om (3.18) 


= . ree... 9: ms | 
3. For anv je L fet v, x’. x, m’ be the sequences and Vs) ania commie. ayemen 
instance Wo= {¥). ¥). .. ¥.J. With this definition the Itkelihood function 


| 


Equation 3.9 associated with the estimate of x satisfies the recursion 





; ae { ; 
+1, 5+ a +iy2 2 fp : 
OT (nm! pe Gh) sep (Se ug nee Lye 4 By(m, 4) CO) 
In the case of the index set is L = {(0,1,....N} the likelihood finetion Eq aime 


3.9 for the estimate x(m) is 
(&%m)) = E(X(m)) (3.20) 


By these considerations the smoothing problem of the ‘data 1s educedmge 
maximization of the Itkelihood function Equation 3.9 with respect to all binary 
sequences m. This can be done using the recursive formula Equation 3.19 and 


Dvynanuc Programming techniques given the details in Appendix A. 


C. SMOOTHING OF PIECEWISE CONSTANT SIGNALS IN TWO 
DINIENSIONS USING KALMAN FILTERING 


There is a problem in extending Kalman Filtering techniques from one dimension 
Into two dimensions due to the lack of a causal state-space model for higher 
dimensions. Anyway, we can determine a two dimensional recursive formulation like 
Equation 3.13 - 3.15 by assuming the estimate X as the average value of the noisv data 
within the regions. 

Proceeding just like the case tn one dimension, a smoothing algorithm can be 


developed to maximize the likelihood function given below: 








N N A 
(k= Pree ¥(k,t)-x(k,t)]? + Bd SX, eX py) + BE Xe Xe) (3.21) 
ket k,t K,t 
To do this, first assume the edges to be known, and define M = ({m, ,} where m,, € 


(€g.€).e5.€,} for each (k.t) € L. e, being the four combinations of edges corresponding 


to the four possibilities in clique n, ' as 


ee) — =! ees) 
and 
ee) 5) 


These combinations of edges is shown in Figure 3.2. 
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Figure 3.2 Illustration of Edges. 


Some definitions are given below referring to Figure 3.3 for each point (k,t) € L to 
obtain the smoothing algorithm: 

Definition 1: t,, 18 the number of points between (Kt) and the closest edge 
above. 


ed): A : , : : : 
Definition 2: Z, , 1s the average value of the noisy data vy in the (t, , * 1) region. 


ht 
Definition 3: he is the number of points i the homogeneous region surrounded 
bv row kK. columin t and the line of edges. 
. e a = . . ° . . ' a . 
Definition 4: X, , 1s the average of noisy data y in the region in which A, | 1s. 


Using the definitions given above, the filtering algorithm can be developed as 


following 
ne ay _ OA oe 
“Kt CK edt | Re AS Seer) (ot) 








Figure 3.3 Illustration of Parameters of Two Dimensional Kalman Fitter. 


where 
? l ae 
Se oy es (3.258 
Kat 
and 
~ a a LH “~ nN ee 
<a et AVA NS era) (3.20) 
where 
t 
as heal ae 
ar t (372m 
k,t 
also 


f 1 
t = ey an ifm. € G : Sie ae 


Kt | otherwise 


a ‘sy eat Sie 


Th otherwise 


rf is 
if My, Ee 


eee 


This algorithm ts based on the fact that for anv disjoint pair of regions A,B © L, the 
Meenas Ol a 101s sicnal denoted by x on the regions :.\.B.,A VU B are related as 

a. “A A — 

mA Ul B= a, X(A) + a,x(B) (3.30) 


with Aes | bl, aa eA UB) syhere 





.| denotes the number of elements in 


fe) 


the set. The recursion Equation 3.24 computes the average on A and Equation 3.26 
tieeaverage On :\ U B. Asmithe One Gimensional case Equations 3.21 - 3.24 represent 
a Well-defined mapping which associates an estimate UML) to any field of edges M é 
Sasa 


(Ege 5.09} Analogous considerations as for the one dimensional case lead 


to the likelihood function as 


A ox ; A ; ae 
Cix{miv, ft = 0% S2(XLn(v, 1} ee 
Waiere 
+ | i me 
mV, ) O if Mer = ey 1 ey (ope) 
-| Vy =e 


So, the estimate x can be determined by searching over all posstble sets of edges M in 
“A , ee ; ; A 

order to find x which maximuzes the Irkelihood function Equation 3.21. However the 

Dvynanuc Programnung algorithm 1s used to maxinuze the likelrhood function in two 


dimensions. 


a) 
Car 
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IV; RESULTS OF SIMELATIOCNS 


A. ONE DIMENSIONAL SIGNALS 

We started the sunulations bv creating a one dimensional two level binarv 
sequence which has 128 points as shown in Figure 4.1. Gaussian zero mean noise 
with different signal to noise ratios (SNRs) catculated” with tie formula en ene 
Equation 3.2 was generated and added to the original signal bv using [NISL horam 
subroutines. The noisy sequences are given in Figure 4.2, 4.5, 4.4 and 4.5. Here. the 


standard deviations for the noise are 25. 50, 75 and 100 respectively. 





—————— ee a mae a -eas aee a ee EE EEE = 


rigure 4.1 Original Binary Sequence. 





Figure 4.2 Noisv Sequence with o= 25 
SNR = 4. 
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Figure 4.3 Noisy Sequence with o= 50 
SNC 2. 
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Figure 4.4 Noiss Seyttence with o= 75 
SV Anta aioe 
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Figure 4.5 Noisy Sequence with ¢ = 100 
SNR = 1. 

The segmentation algorithm explained in the first section of Chapter 3 which uses 
Dvnamic Programming was applied to filter out the noisy signal with the signal to 
noise ratios given above. To tmprove the result, for each simulation the agormhm 
Was run twice, from left to right and from right to left in order to provide sinooth 
results. 

The outcome of simulations depend on the model parameter J that models the 
smoothness of the original data and this parameter was switched several times until we 
had the best segmentation for each case. The values of Pf were Set Dv trial ade uen 
for this study. The best values for B for the given SNRs were obtained as 0.34. 9.80, 
O.74 and 0.72 respectively. Filtered signals are presented in Figures 4.0. 4.7, 4.8 and 


eo 


30) 


Figure 4.6 Segmented Sequence for o= 25 and p 


O.S4. 


Figure 4.7 
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Segmented Sequence for 0 = 30 and B = V.S0. 
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Figure 4.9 Segmented Sequence for ¢= 100 and B = U.72. 


for the case of data with unknown levels we developed a different secimentagienm 
algorithm. The difficultv 1s to estinate the dara levels, as well as to determinewale 
segmentation. This segmentation algorithm is based on Kalman Filtering and 
Dvnamuc Programming. First. this method has been applied to the one dimensional 
noisy signals shown in Figures 4.2, 4.3. 4.4 and 4.5. The noisy data was filtered 
twice (two passes). The two passes results are shown in Figures 4.10, 4.11 for o= 25 


and 30, 


Figcure 4.10 Segmented Sequence by AF with a= 25. 
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Figure 4.11 Segmented Seauence by RF with o =v. 


B. TWO DIMENSIONAL IMAGES 

We applied the two dimensional filtering algorithm to the test images which 
contain three different resions given in Figure 12>" im this tesmimage thicre anemic 
different intensitv levels (one for the background and three for the objects). Random 
Gaussian noise has been added with standard deviation o= 10 and o= 20 as shown in 
Figures 4.13, 4.14. 

The result of filtering is shown in Figures 4.15 and 4.16. Notice that in these 
cases the noise has been removed while preserving the edges between the regions. The 
values of the parameter B giving the best results are B=2.0 for o= 10 and B=1.S for 
6= 20 showing a trend of this algorithm. The value of B depends on the SNR of the 


noisy picture. In particular } should be decreased for hugher levels of noise. 





Sea mecppicaion (ordetest imace Im Wnderwater environment, We apolied mie 
meen to @ 512 < 512 picture of a fish) corrupted by additive noise. The nos 
Image is shown in Figure 4.17 and the filtered one is shown in Figure 4.18. The 
significance of it is the fact that after filtering the fish and the background present weil 
distinet intensity levels. Although applications to tus class of problems are sull under 
Misra, tlie fact hat the Object (the fisli im this case) presents characterisucs well 
distinet from the background can be used for automatic detection er recogniuion m an 


artilicial intelligence framework. 





Figure 4.12 Ortginal Test Image. 





Figure 4.13 Test Image with o= 10. 





Figure 4.14 Test Image with o = 2u. 
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Figure 4.15 


Segmented [mage with o= 10 and B= 2.0. 
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Figure 4.16 Segmented Test Image with o= 20 and B=LS. 
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Figure 4.18 Segmented Fish with o = 25 and B=0.235. 


Analogously we tested the algorithm on a 16 level checkboard shown in Figure 
4.19 for different levels of noise. In particular the values of the 16 levels are given by 
Ne lee Oe OU 120, 60, 10, 75, 175, 90, Gd, 130, 55, 190, 110 from ictt to right 
for each line. The effects of additive noise are given in Figures 4.20 and 4.21 tor o= 10 
and o = 20 respectively. The noise 1s alwavs assumed to be Gaussian. zero mean ind 
White. The application of the filtering algorithm is shown in Figures 4.22 and 4.25 
using B=1.2 and B=0.7 respectively. As expected the algorithm filters within the 
regions While preserving the sharp separation between adjacent regions. .Also shown in 
Figure 4.22 are the edges between regions, as detected by the algorithm which shows 


the reinitializations of the Kalman Filter gains. 





Figure 4.19 Test Image with 16 Ditterent Regions. 
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Pigtines 20esiliest Image with o= 10. 
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Figure 4.21 Test Image with o= 20. 
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Figure 4.23 Segmented Image with o= 20 and Bp =0.7. 
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y: GONCEUSIONS 


An algorithin to segment noisy data in one dimension and two dimensions hus 
Beewepre-mred, Wiicice the Seements are piecewise Constant. Ihe data are described bv 
state-space models. The particular feature of the algorithin is to search for the best 
sequence of edges in order to mantmze the lkehhood function. The algorthm. has 
emphasized piecewise constant data, but it can be extended for the data with iegions 
characterized by textures to Which we associate different Autoregressive (AR) models. 
Typical applications are not only image segmentation but also segmentation in speech 


analvsis, such as phenomenon separation. 
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APPENDIN 3 
DYNAMIC PROGRAMMING ALGORITHM 


The Dynamic Programmung ts an alygorithni used when the solution to a problem 
mav be considered as the result of a sequence of decisions. An optunal sequence of 
decisions will maximize the given function, in the case that 1s used in this study, the 
likelihood function €,. 

In Dynamic Programming an optimal seyuence of decisions 1s obtained bv 
making explicit appeal to the Principle of Optimality. In the cases when this principle 
can be applied. an optimal seyuence of decisions has the property that whaiever the 
Injual state and decision are, the rest of the decisions must make up an optimal 
decision sequence according to the state resulting from the first decision. Standard 
Dynamic Programming techniques are introduced by Bellman. [Ref. 9] 

The mathematical equations and the algonthin is given below for the case that 
was mentioned in the second section of Chapter 3. By defining the likeithood function 


ats) 





k ik 
= ? 
E(XgsX prey) = BY 865.811) - hy FOI (A.1) 
j}=0 207 i=0 
Or in recursive form 
er Bry cre thse y= 2G ecg J PAN Gen oo TA) (A.2) 


where 





7 Pe 
[Scere ene Reus.) 


2 


AE OX eb 1X) = BEC 4 9X) - : 


and setting the initial condition x_, = 0 and also assuming that x, can be only logic 
"O° or 1 we can, determine theroese seduce: iy , kK=0,....N} Which maxinuzes 


E\(Xy,.4XX). For this purpose define 


J(O) = max SG: ar Xe oy.) ee) 
oN? 
and 
eax Bde Xo) Nes 


(SgeesX 1) 


Weeweat) represent the proolem in the form of a graph as in Tigure A.t. with N stages 
Where the nodes at the k-th stage represent the state x, (either “O” or “17i, and the 
eenes Fepresent the wal associated to the tramsition from one state to the next. In 


this way J,(0) and J.(1) represent the maximum of alt possible gains up to x, = U and 


x, = Ll respectively. These quantities can be recursively updated as 

i ey) = max; de (0) + AE 4 mee Jie Nee (1,09) (A\.6) 
and 

Poa = maxis (0) + a ee 00 UE er yer aN feel (A.7) 


with AC, as in Equation A.3. Also we can keep track of the branches which vield the 
maximum likelihood by defining Pointer,(0), and Pointer,(1) at each stage kK. The 
sequence X maxinizing UN ii tnw) is therefore obtained bv backtracking as in the 
following: 

eet X be such that INR.) = maxJ\(0) , Jx( Le 

for kK. — \ =ihte 0 
Ne = Rombcreeants,. ) 
end for 


end 


Ol 





k=0 


logic "0° 


logic "1° 


Figure A.1 


Choice of Best Sequence of Decisions. 
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APPENDIN B 
COMPUTER PROGRAMNIS 


KRREEKARRERAKRKEKARKAKKRKRRARKRKKRKKRARRKRKRRRKRKRKRKARKARKRKRKKKRKRRRKKKRRRARRARE 


PROGRAM BINARY.DAT 

PURPOSE To generate 128 point binary sequence. 
REQUIRED IMSL ROUTINES None. 

DIESE iets, Lerehdild HAclPAsAOGLU Jan.71987 


ee Ne ee GC ee RT RRR RR IRE KK RARE RR AR RAKE 


AA A A A OO eA 


x 
x 
* 
x 
* 
* 
x 
x 
* 
* 


Pigteeee Tiltat FINI P(i26) TERI, SETBUP ,K 

REAL DWINDO,MNOVEA, DRAWA 

OPEN(UNIT=1 ,NANE='BINARY.DAT' ,ACCESS='DIRECT' ,STATUS='NEW' , 
REGE=32 VARREC=12e) 


DO 10 K=1 026 
Tr GK iia Plone aid. 30') THEN 


Fick 

Pese ew sclamy Or .AND. Kk .CT. 90) THEN 
F(K)=100 

Poe wnGlent ZO AND. K JLEe i266) THEN 
F(K)=100 


oe 
F(K)=200 


END IF 
WRITE(1'K) (F(K)) 
CONTINUE 

CLOSE(UNIT=1 ) 


CALL INITT(480) 

CALL TERM(2,1 

CAEL Sag 

CALL DWINDO(-50.0, 300. 50 ,750-0,400.0) 
CALL MOVEA(-50.0,-50 

CALL DRAWA(300.0,-50. i 

CALL DRAWA(300.0,400.0 

CALL DRAWA(-50.0,400.0) 

CALL DRAWA(-50.0,-50.0) 

CALL MOVEA(0.0,1000.0) 


CALL DRAWA(0.0,0.0) 

CALL DRAWA(128.0,0.0) 

CALL MOVEA(0.0,0.0) 

DO 35 K=1,128 

K=K 
=F (K) 

CALL DRAWA(X,Y) 
CONTINUE 

CALL FINITT(0,0) 

STOP 


END 


AAO AAAANOAOMM1© 


iC 
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KRAKKKAKAKAKKKKAKKAKKKKAKAKKARAKAKKAKAKKRKAKKKKAKRAREKRRERE 


PROGRAM BINARY2. DAT 


PURPOSE To generate one dimensional 128 point binary 


sequence 


with zero mean Gaussian noise used 


Cue Ls ys Leila 


REQUIRED IMSL ROUT INES  Geliiie 


*~ 

* 

* 

~ 

* 

x by Dynamic ee ens algorithm to segment 
x 

* 

* 

* INPLEMENTED BY Lt. 
x 

~ 


Kani HACIPASAOGLU Feb. 1987 


INTEGER INZTIT FIN(@@, F( 128) teh Sere wae eee 
REAL DWIHDO, MOVEA, DRAWA, GGNHL , R, $(128) 


DOUBLE PRECISION DSEED 
DSEED=65471 


OPEN (UNIT=2 ,NANE='BINARY2.DAT' ,ACCESS='DIRECT' ,STATUS='NEW', 


* RECL=32,NAXREC=128) 


DOmwLO I=1 7 tzs 
PEAT GE. Ome bee 
=109 


Lie sO )miaeh 


rey Ba oe 70 .AND. I .LT. 90) THEN 
ELSE IF(i Gre 120 SAND, Geeke los memes 
F(I)=100 


ELSE 
F(I)=200 
END IF 
NR=1 
CALE GGMMIL (DSEED, NR,R) 
SIGNA=75. 
P= eTCHASR 


Pe NOUSE eae WHICH WAS GENERATED 


S(I)=F(I)+R 
URTIB (20h) (sia 
10 CONTINUE 
CLOSE (UNIT=2) 


CALL INITT(480) 
CALE SE RuteZ ra: 
CALL Set eUe 


2 
CALL PE Ngo Gan) 200.3 


CALL HOVEA=S070> = 5020 
CALL DRAWA( 300707 5070 
CALL DRAWA(300.0,400.0 
CALL DRAWA(-50.0,400.0 


CALL DRAWAl 5070-5020 


CALL MOVEA(0.0,1000.0) 
CALL DRAWA(0.0,0.0 
CALL DRAWA(128.0,0.0) 
CALL MOVEA(0.0,0.0) 
DO 35 I=l, 128 

X=I 

Y=5(T) 
CALL DRAWA(X,Y) 

35 CONTINUE 

CALE FINE UE Coa) 
STOP 
END 


,- 50. 0540020); 


Cs 
AES 


20% Pas doo ES SE OS 


' ’ 


OAD OR 


MOAANDNDANUAANANAN 


FROGRAM ONUR2 


PURE OSE (eye -cuienibea olay binaky sequence using 


x 
x 
ZS 
- 
Dynamic Programming algorithm. 
: REQGULTRED IiSi ROUTINES iene. 

x 


IMPLENENTED BY Lt. Kani HACIPASAOGLU March 1967 


PUSECER eke bOnr] OULsPIR( 128,03 4), F(128) 
MIEECER SE FINITT, TERI, SETBEUF, | 


REAL BETA,SIGiIa, LHIENO, LNEWL,NOVEA, ‘to. Pigeon Prt2oy sult, TsE , 10 
REAL DELV(1:128,0:1,0:1),DWINDO, DRAWA,S(128) ,X,Y,SSE,SEE, EE 


SIGHA=159.0 
FAD(*,1) BETA 

roman (Fa. i» 

FO=10 

PLS 300 
OPEN(UNIT=2 ,NAME='BINARY2.DAT!,ACCESS='DIRECT', 

* RECL=32,MAXPEC=128) 

Poros ti=1 ize 
READ (2'N) (S(N)) 
WRITE(6,%*) S(N),N 


35 COL BiUE 


rs 


CLOSE (UNIT=2) 
K=1 


SHatus= OLD". 


SiGiA. 22 
SLeMae 2) 


(SIGIA**2)) 
( 
/1,0))) THEN 


) 
(stain) 
h 


ee BO OG ee) 
L1=-BETA-(((S$(K)-200)**2)/(2*(SIGHA**2) )) 
DO 15 N=1,128 
K=M-1 
DELV(K+1,0,0)=BETA-(((S(K+1)-FO)**2) ee 
DELV(K+1,1,0)=-BETA- rt} -FO ee ea 
DELV(K+1,0,1)=-BETA- Red -Fl ee Ke 
DELV(K+1,1,1)=BETA-( S(K+1) 5 F1)**2)/(2% 
IF ((LO+DELV(K+1,0,0 (L1+DELV(K+1 
LNEWO=L1+DELV(K+1,1 oy 
PTR(M,O)=1 
LSE 
LNEWO=LO+DELV(K+1,0,0) 
PTR(IM,0)=0 
END IF 
IF((LO+DELV(K+1,0,1)) .LT. (L1+DELV(K+1,1,1))) THEN 
LNEW1=L1+DELV(K+1,1,1) 
PTR(M,1)=1 
ELSE 
LNEW1= LO+DELV(K+1, eels) 
PTR(M,1)=0 
END IF 
LO=LNEWO 
L1=LNEW1 
CONTINUE 


TRACE BACK PROCEDURE 


Clee 3  NeME= OUIS.,DAT*,ACCESS='DIRECT' ,STATUS='NEW', 


* RECL=32,MAXREC=128) 
ImeO> Li. Li) THEN 
J=1 


Oug=J 
DO 44 K=127,0,-1 
IF(OUT .EQ. 0) THEN 
OUTS (K+i )=100 
ELSE 
OUTS (K+1)=200 
END IF 


ta 
Ca 


OUT=PTR(K+1 , OUT) 
WRITE(6,*) OUTS(K+1),K+l 
44 CONTINUE 
DO 55 K=1,128 
WRITE(3'K)(OUTS(K) ) 
55 CONTINUE 
CLOSE (UNIT=3) 


CALL INITT(480) 

CALE TERING! 

CALEP Ss ee 

CALL DVWIN Su 390. OF =50 707 400. Up 
CALE HOVE] S070) - 0-0 

CALL DRAWAV SCO 0 72070 

CALL DRAWA 300.0,400.0 

CALL DRAWA(-50.0,400.0 

CATE RRR eGo. 0) SoU. 


0. 
CALL DRAWA(150.0,0.0) 
CALL MOVEA(0.0,0.0) 


DO 3 K=1,128 
ik 
Y=OUTS (K) 
CALL DRAWA(X,Y) 
3 CONTINUE 
CALL FENZTECO 0G) 
STOP 
END 


AAA GAHAAAAAANAON 


Pte ai, WP a, Ga, ln fo, Ga, ame, Gr | 


RAKKRKAKKKKKAARAKRKRKRKAKKKAKKKAKKKRAKKKRARRERKRKKARRARARARM AERA 


PROGRAM TEST.DAT 


PURPOSE @iicmceneraccmantesc Image whach has 
three different intensity levels. 


REQUIRED TSE ROUBEVE Sa tone . 
TNPCE MENTE er Ue rAeCIPASAOQCLU) May 1967 


Bee A(lZ6,l2e) 


x 
* 
~< 
~ 
x 
* 
* 
as 
x 
i 
& 


wee Chr shy AG lmniele, Ke@2zpcec na lel Ze: 128 )eb 2 leks, 12123 ) 


R=20 

AC1=60 
YC1=60 
YC2=60 
XC2=90 


OPEN(UNI1=5 MeMe-"LEST.BAT' ,ACGESS='DIRECT' ,STATUS='NEW' , 


RECL=32 ,MAXREC=128) 


DO 10 I=1,128 
DOYZO07 9-1 Zo 
Fl Esha UE HCa) 25 (3-22) 42) - (Rend ) 


Pole) T-XC2 )**2+( J-YC2 )**2 R**2)) 
ua eM SL Rroeeopmeano. (f2(1.3) .Gr. 0)) THEN 
ACL. we 
ance a) Fit, J) Rem Oy ee AND. Ch20Tu) LE. 0)) THEN 
I,J)=1 
ELSE IF (FL(T, J) Reo ye AND. (rPecneo)) .LE. 0)) THEN 
j 


ELS 
A(I,J)=200 
END If 


GOMEINUE 


WRITE(5"I) (A(I,J3),J21,128) 
UE 


CONTIN 


CLOSE(UNIT=5) 
=) Oks 
0 b8, 
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PROGRAM  TESNENSDAT 


- 
- 

x 

x PURFOSE To generate a test image with zero mean 
z Gaussian noise to be segmented by 

i. Kalman Filter. 
* 
* 
* 


REQUIRED SIMSEMROURENE Ss GGUIE 
IMPLEMENTED BY Lt. Kani HACIPASAOGLU May 1987 


* 
nx 
RAKRAKKKKKRARKRKARKKRR RE KKK RARKRKRKRKRRRKKRKRRRR RRR RRKKRKKRRARAR RARER 
BYTES AZo ce) 
CHARACTER*50 EES i eh 
INTECER@ Grol) necro, F1( 128 428) 
INDEG ER NR, Rr, SIGMA, Ck F2(128 128) ,AA(12Z3 ize 
REAL GGNUL,R 
DOUBLE PRECISION DSEED 
DSEED=65471 
WRLIEt A700) 
100 FORMAT('ENTER IMAGE DATA NAME' ) 
READ Gol 
101 FORMAT CAS©) 


CR=20 

KC1=60 
YC1=60 
YT2=60 
XC2=90 


OPEN(UNIT=1 ,NAME='TESTIM.DAT' ,ACCESS='DIRECT' ,STATUS='NEw', 
* RECL=32 ,MAXREC=128) 


WRITE(*, 222) 
222 FORHAT ( | ENTER SIGMA!) 

READ(*,333) SIGMA 
333 FORMAT (12) 

DO 10 [= 22 

DO 20 J=1,128 

ES ee oe SEs yee Al CR**2 
E2(h I-KC2 )**24+( J-YC2)**2 ae 
TF ((FL(Z J) ~LE. ©) .ANDs (CF2( 130) .Gl.. ©) mine 


— TEC(EL(T, J) -.GT.9O) J2ND. (h2(1sd DEO pe rhe 

A(I 

ELSE CLC, 3) EE. 0) AND. (F217) Lene mea 
Dies 

S 


LSE 
A(I,J)=200 
END IF 
NR=1 
CALL GGNML(DSEED 
Rr=SIGMA*jnint(R) 
AAT “J)sAlL, J)+Rr 
AA(I,J) .LT. -128) THEN 


AA ( \= +256 
ELSE BENG cp GT 327) THEN 
AA(I J)=AA(I, nes 


Py see 

END IF 
20 CONTINUE 

WRITE(1'1L) (A(1I,J),J=1,128) 
10 CONTINUE 

CLOSE(UNIT=1) 

SfOr 

END 


NR,R) 


OY OAOAADAAAAAAGG 


20 
10 


30 


* 
x 
ae 
* 
x 
x 
* 
“< 
= 
a 
* 


PROGRAM 
PURPOSE sae 


qin 
REQUIRED IMSL ROUT Iies ftlone. 
Pie LEMENTED Sa aa is. 


Beto ACiZe 023) 


Cee UNM=lalip= “[iALsoat!  ACGESS='DIRECT" 


i 


FINAL .DAT 


enerate a test image which has ls 


erent regions. 


ReGL=3 2 lich eo= 126) 


DO 10. .K=1,32 
BO 20) Lele oz 











ae sae 
Miku oe2 =soO 
eee 
A(K,L+96 )=70 
ihe 200 
A(K+32,L+32)=120 
A(K+32,L+64)=60 
A(K+32,L+96)=140 
A(K+64,L)=75 
A(K+64,L+32)=175 
At K+64 , L+64)=99 
A(K+64,L+96)=65 
A(K+96,L)=130 
A(K+96,L+32)=55 
ies =190 
A(K+96,L+96)=110 
CONTINUE 
CONTINUE 
Homes 1-112 


RITE(L"T) (Ale) = 123) 
COMLINUE 


CLOSE(UNIT=1) 
STOP 
END 


Kanl HACIPASAOGLU June 1987 


eae ie ae eae ee et 


Cie 2 Ue ace eee 


pole S—" HEN", 


AAMAANANANANAANN 


[Tee Geen ee RS aro 


To caeuace a test image which has 14 
dia 


erent regions with zero mean Gaussian 


BYTE SAC a2erizs) 


* 

- PROGRAM NFINAL.DAT 
- PURPOSE 

Pad 

* 

a REQUIRED 

x 

cs IMPLEMENTED BY Lt. 
x 


IMSL ROUTINES GGNHML 
Kani HACIPASAOGLU June 1987 


a 
x 
noise to be segmented by Kalman Filter. a 
« 
* 
Fas 





CHARACTER*59 NFINAL 
INTEGER NR,Rr,SIGMA,AA(128,128) 
REAL GGNHL,R 
DOURGE@E Roel a LONe Porc) 
DSEED=65471 
leh sels (es LOO 
100 FORMAT( ' ENTER IMAGE DATA NAME' ) 
ReaD’ Lo) 
101 FORMAT(A50) 
OPEN(UNIT=1 ,NANE='NFINAL3.DAT' ,ACCESS='DIRECT' ,STATUS='NEW' , 
A RECL=32 ,MAXREC=128) 
WRITE(*,222) 
222 FORMNAT( ' ENTER SIGMA!) 
READ(*,333) SIGMA 
333 FORMAT(I2) 
DOe1O Kal) 32 
DO 205 wale 32 
A(K,L)=1900 
A(K,L+32)=590 
A(K,L+64)=180 
A(K, L+96)=70 
A(K+32,L)=200 
A(K+32,L+32)=120 
A(K+32,L+64 )=60 
A(K+32,L+96 )=140 
A(K+64,L)=75 
A(K+64,L+32)=175 
a Kred +4 =90 
A(K+64,L+96 )=65 
A(K+96,L)=130 
A(K+96 ,L+32)=55 
A(K+96,L+64)=190 
A(K+96,L+96 )=110 
20 CONTINUE 
10 CONTINUE 
NR=1 
DO 80 I=1,128 
DO 90 J=1,128 
CALL GGNML(DSEED,NR,R) 
Rr=SIGMA*jnint(R) 
BA(T J)=a(1 J) +Re 
TP Ga( ie) eee — Zo eee 
(I,J)=AA(1I,J)+256 
ELSE IF(AA(I,J) .GT. 127) THEN 
I,J)=AA(I,J)-256 
Cie =n ae) 
END IF 
A(I,J)=AA(I,J) 
30 CONTINUE 


WRITE(1'I) (A(L,J),J=1,128) 


60 


80 


CONTINUE 
SpooECUNIT=1) 
SOP 


END 


61 
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PROGRAM “GR iol reD 


PURPOSE To seament a two dimensional noisy image 
using Markov Random Field, Pynamic 


REQULTRED ST SheROUTINESs leme. 
THPLEMENTED BY Lt. Kani HACIPASAOCEU MeAticnm 57 


x 
x 
” 
x 
~ 
; Programming and Kalman Filtering techniques. 
x 
x 
x 
x 
x 


RAEEARAKRKKKKRKAKRKRKRKRKRKRRKARKRKRRERRKRKRRKERARKRRKARRRRRRKKERRRRRARERRARER? 


commoneny, visolz 512), tau(sl2) )xtop( size tedge, ac. 
integer start, dir, cycle 


character*20 image, imout 
byte Ge ee ioe -ol2) 


write ( 

format (' Enter Input Image File Name:') 
read (*, image 
write(*, 782 

format(! Enter Output Image File Name:') 
read(*, 778) imout 

format(a20) 


open (2, name=image, access='direct', status='old') 
get data from file 

write (*, 779 
format | Enter Image Size') 

read(*, 780) npoints 
format (i3) 

Be). oe . . 

Sloe ee eal 
write (*, 790 - 


format(i1) 
oe ; 20), 
write(*, 792 
pr Enter scale factor (1.0=no scale) ' ) 
format output yout= =yth+scale*(y- Vieni 
read(*, 793) scale 
formal floes 
write ( hie 
format (| pone ian) 2 2:yth=0.9ymax; 3: enter yth; 
read(* 1yth 
lf 7 ee eq.3) then 
mrss. 7 4) 


formatt eenber yen» 
paca (23) yth 


do 5 k=l, npoint | . 
read(2") (S (ie, j), j=1, npoints) 
do j=1, my eh Ss 
temp Pie, d? 
ah eonee t.0,0) then 
temp=temp+256 
endif 
y(k,j)=temp 
continue 
continue 
close (unit=2) 


pena ENTER: sigma, beta') 
format (210.4) 


asters) jiilel Liam Korey ofa tpicass 


? . ' 


x 
* 
x 
ree 
~ 
x 
2 
ra 
4 
x 
Pas 


fare) 


Aan 


write(*, 667) 
foumat(' StarceGonpucing sae) 
ymax=9.0 
yav=0.0 
ncount=0 
nx=noolnts 
ny=npoints 
gam=1.0/(2.0*sigma**2) 


Roe t Lal eae 


do 30 j=l,ny 
tau(})=0.0 
ease Oo 
Geneinue 
nX=nx-1 


eect Loop  —~s> 
do 40 1=1,nx 
Stant=| 
dir=+1 
cycle=l 
call line(i, start, dir, cycle, gam, beta, rowav) 


start=ny 
cycle=2 
call line(i,start, dir, cycle, gam, beta,rowav) 
COMmt Live 
x** compute average of output image 
yav=(ncount*yav+rowav)/(ncounttl 
neount=neounttl 


write(*, 668) 


EGmiat(’ “Wmeonad Commuting. NOwmeranstering data to disk ... 


ARKX output data *** 


if (1ytheege2)) then 
yth=ymax*0 .9 
. eeenait 
if (A iyehseq 2b) tien 
yth=(yavt+ymax)/2.0 
endif 


UrLite =, 800) 
format(' yave, ymax, yth =') 
write(*, 801) yav, ymax, yth 
format(3(x, £10.4)) 


do 60 k=l, npoints 
Gomel j=l, tpoints . 
ys=yth+scale*(y(k,j)-yth) 
ies Guy eos.0) then 
ss 0 
endif 
tee ys. be. 0.0) then 
s=0.0 
endif 
itemp=jnint(ys) 


ifeereenp.gde.i127) then 
itemp=itemp~256 


endif 
v(k,j)=1temp 
GOntlnue 
continue 


') 


open (unit=3, name=imout, access='direct', 
* Status='new!', recl=128, maxrec=512) 


do 50 ae Npo une { 

write ( Coie on k=1,npoints) 
310 continue 

close =a 


stop 
end 


subroutine alate start, dir, cycle, gam, beta, rowav) 
GIG (aie Species ieec yc. ~~ 

common ny, y(° 25 be. SS ae xtop(512),1edge, ymax 
real ate IZ, Lambda(4,5 

real ttau ae txtop(4), Tee G) txhat(4), te(4) 
real Bea 512)" newxtop (4,5 2 

real e(4,5 

integer Bean a. Bl ) 


e *xAX scan from START by DIR (=+1, -1) 
jEstart 

Cc AAR initialize first column 
do 90 k=1,4 
xhat«k,startj=y7(i,start,) 
lamb a(k, Shaye =0.0 
e(k, start)=0 

90 continue 


rowav=0 .0 
heount=0 


i5jtdir 
ne Statea Sac fee ed eu upper edge) 
xx QO=no e 


kX state (0,0) 
do sO k=, i 
Sas tau( 
Exto Fs )=xtop 3) Te Ci, Sediryeeeautk) ,j)*xtop(j)) 
tlambda(k) Se with j-dir)+ttau 

Z=h eee (oie smcbie, 9 - 
txhat(i)= aCe j- airy + (tau Od /LleRoeeere 


el= sel bee Se ae 
e2= 1,j)-txhat(k) 

iti,j)- cated )) 2 
sat e(k, } ~dir)+gam*(el+e2+e3)/3 - 2%*beta 
if (k.eq.1) then 
emin=te (i) 


fe Geauevelenentin then 
im=k 
emin=te(k) 
ence 
IEG) continue 


gana 


newtau 


lambda )=tlanpea 


=tlambda(im) 


er 


eee (im) 
xhat(4 |) =txtop | 


e(4,j)= 


Axx state (0, fy KKK 
do 115 k=] 
tlambda(k)= Be eaae Bells >=) 
z=y(1,j))-xhat(k,j-dir) 


To 
O 
-- 
— 
ct 
@m 
tae XY 
crite ~_ Bot 
O~— Pp fp~ ™~ 


a0 
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txhat(k)=xhat(k,j-dir)+(1.0/tlambda(k))*z 





@ 
as Yi} pCa eka hig 2 
zaraee Ne tahatliey 2, 
ie j)-txhat{ “)) ae 
oO =e(k,j-dir)+gam*(el+e2+e3)/3 
C 
if (k.eq.1) then 
ee te(l) 
eis olog 
| BECK ie, emin) then 
im=k 
emin=te(k) 
. endif 
115 Gene lune 
fe 
newtau(1,j)=1 
newxtop(1, la 
lambda(1,j ep cee 
Bearer a on 
lsh el Guay) Heed keer obiuly, 
e(1,j))=te(im 
Cc 
c *kA state ex fe 
do 120 k=l, 
Seas Scare ‘ae . 
. Geren = ceo Gg ete s(t .0/ttau(k))=(¥(1,3)-xtop(j)) 
el=(y(1,jtdir)-txto my) 
eae ee sj) stweoo eka 
e3= he: We exuop.) )* 
te ie =e(k,j-dir) + SE eee 
c 
if (k.eq.1) then 
ee 
ie (tee) ere.emin) ten 
im=k 
emin=te(k) 
. enaue 
120 continue 
Cc 
Cc 
newtau(2,j)=ttau(im) 
newxtop(2 *4)= txtop(im) 
lambda(2,j)=tau G) 
seemcon (2 4) =im 
xhat(2,))=txtop(im) 
e(2,))=te(im 
C 
@ 
C Eeestate <( 1), 1) 
ee (44) 
NeW ceOpVo,))/-Y\1,) 
lambda(3 ven . 
xhat (3,3) y( a 
Sc gee 
es=(y1t1,3)- S534 az 
z2= ee lelte3) S 
moe Jeet). toe i jedir)) then 
e(3, j)= z2+e(4, j-dir)+beta 
Boniiee nic) 3} )=4 
else 
e(3,j))=z2te(l1,j-dir)+beta 
. pointer (3, j)=1 
endif 
Cc 
C 


Memory (l-dat).and.j.1t.(my-dir)) then 


65 


aAaaNAA 


an a0 


an 


GO un 72 
endif 


AAK DaCwerdG anaes 


emin=e(1,}) 
do 150 k=1,4 | 
if (emin.le.e(k,j)) then 
Hanes (ae 
im=k 
. endif 
continue 


y(i,j)=xhat(im,j) 


WP te (4 alee) aca 
format (is) 


if (cycle.eq.2) then — 
xtop(j)=newxtop(im,j) 
SER cuca Cimiast 


APASEGMIPUCe Max wiMeetio dey 
if (y(i,j).gt.ymax) then 


— ymax=y(1,3) 
ee 


xx row average =. 
rowav=(ncount*rowavty(i,j))/(ncountt1l) 
ncount=ncounttl 


mark the edge with a black entry 
tf Mmene 4 cid; tedde cece amemon 
Gy) 00 
enol is 
ailolahag 


im=pointer(im,j) 

=) cnn 

i¢ (GaGgt.(l=dir).end.).1t. (ny-darwmenen 
goto 888 

endif 


return 
end 
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